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With the goal of taking a step toward the construction of astrophysically realistic initial data for 
numerical simulations of black holes, we for the first time derive a family of fully general relativistic 
initial data based on post-2-Newtonian expansions of the 3-metric and extrinsic curvature without 
spin. It is expected that such initial data provide a direct connection with the early inspiral phase 
of the binary system. We discuss a straightforward numerical implementation, which is based on a 
generalized puncture method. Furthermore, we suggest a method to address some of the inherent 
ambiguity in mapping post-Newtonian data onto a solution of the general relativistic constraints. 
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I. INTRODUCTION 

One of the most exciting scientific objectives of gravi- 
tational wave astronomy involves the search for and de- 
tailed study of signals from sources that contain binary 
black holes. Mergers of two black holes both with masses 
of ~ 10 — lOOAf0 will be observable by the ground based 
gravitational wave detectors, such as GEO600, LIGO and 
others Q . These systems are highly relativistic once they 
enter the sensitive frequency band (~ 50 — 200 Hz) of 
the detector. For LISA, gravitational waves from super- 
massive binary black hole mergers (e.g. black holes with 
mass greater than IO^Mq) are very strong, with high 
signal-to- noise ratios up to 10^ ||^, making these events 
observable from almost anywhere in the universe. As- 
trophysically realistic models of binary black hole coales- 
cence are therefore required to study these phenomena in 
detail ||. 

To solve the full Einstein equations in the dynamic, 
non-linear phase at the end of the binary black hole inspi- 
ral we turn to numerical relativity. Numerical relativity 
has advanced to the point where a time interval of up to 
40M (where M is the total mass) of the merger phase of 
two black holes can be computed if the black holes start 
out close to each other ^ ^. Recent simulations of 
head-on collisions of black holes last significantly longer 
and give reason for optimism for the orbiting case 0| . An 
approach to produce at least moderately accurate mod- 
els for the wave forms generated in binary black hole 
mergers was recently developed in the so-called Lazarus 
project ^, |ll|, a technique that bridges 'close' 
and 'far' limit approximations with full numerical rela- 
tivity. This approach has lead to the first approximate 
theoretical estimates for the gravitational radiation wave 
forms and energy to be expected from the plunge of orbit- 
ing non-spinning binary black holes to coalescence ^, ^ . 

Due to theoretical and numerical limitations, all cur- 
rent numerical simulations must begin by specifying ini- 
tial data when the black holes are already very close (sep- 
aration < 7M). There is a push to place the starting 
point of these simulations at earlier times, say at a few 



orbits before a fiducial innermost stable circular orbit 
(ISCO) which approximately marks the transition from 
the inspiral phase to the plunge and merger. But what- 
ever the starting point, the simulation will only be as- 
trophysically meaningful if it starts with astrophysically 
realistic initial data. 

The question we want to address in this paper is there- 
fore how to obtain astrophysically realistic initial data for 
numerical simulations of binary black hole systems. In 
general relativity the initial data must fulfill constraint 
equations, so only part of the data are freely specifiable, 
and the rest is determined by solving the constraint equa- 
tions (for a review see e.g. |l^). A lot of the work in 
constructing initial data has focused on approaches that 
pick the freely specifiable part of the data with the aim 
of simplifying the constraint equations, rather than us- 
ing astrophysically realistic initial data. A standard as- 
sumption is that the 3-metric is conformally flat and the 
extrinsic curvature is derived from a purely longitudinal 
ansatz (see e.g. [|3[ |l^, 16|). Currently, there are a 
number of new approaches |17, |l^, |l| to specify 

'improved', including non-conformally flat, initial data 
for binary black holes. 

However, none of these approaches to construct ini- 
tial data makes explicit use of information from an ap- 
proximation procedure such as the post-Newtonian (PN) 
method, which is believed to accurately represent astro- 
physical systems in the limit of slow-moving/far-apart 
black holes. An approximate binary black hole metric 
based on post-l-Newtonian (IPN) information in a coro- 
tating gauge has been derived by Alvi [Q. However, at 
present this metric cannot be used in numerical simula- 
tions due to the presence of discontinuities in the match- 
ing regions . An interesting approach based on quasi- 
equilibrium sequences of initial data has been studied nu- 
merically, e.g. |2^, although some aspects of the method 
appear to be based on Newtonian or IPN assumptions. 

In this paper we describe a method to generate new 
fully general relativistic initial data for two inspiraling 
black holes from PN expressions. The motivation for 
this method is that even though PN theory may not be 
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able to evolve two black holes when they get close, it 
can still provide initial data for fully nonlinear numeri- 
cal simulations when we start at a separation where PN 
theory is valid. In particular, we obtain an explicit far 
limit interface for the Lazarus approach. Our method 
allows us to incorporate information from the PN treat- 
ment and should eventually provide a direct connection 
to the inspiral radiation. 

Like in other approaches, we start from expressions 
for the 3-mctric and extrinsic curvature in a convenient 
gauge. We use expressions for the 3-metric and its con- 
jugate momentum up to PN order (w/c)^, computed 
in the canonical formalism of ADM by Jaranowski and 
Schafer ||2^. This order corresponds to 2.5PN in the 3- 
metric and 2PN in the conjugate momentum, since the 
latter contains a time derivative. Therefore, the PN data 
are accurate to 2PN. 

The 3-metric and its conjugate momentum are derived 
together with a two-body Hamiltonian using coordinate 
conditions |2^, which correspond to the ADM 
transverse-traceless (ADMTT) gauge. Note that there 
are several other formulations and gauges for PN theory, 
see e.g. for a review. The ADMTT gauge has several 
advantages: (i) we can easily find expressions for 3-metric 
and extrinsic curvature, (ii) unlike in the harmonic gauge 
no logarithmic divergences appear, (iii) for a single black 
hole the data simply reduce to Schwarzschild in standard 
isotropic coordinates, (iv) up to (v/c)"^ the data look like 
in the puncture approach u&, which simplifies calcula- 
tions, and (v) the trace of the extrinsic curvature van- 
ishes up to order (v/c)^, so that we can set it to zero 
(if we go only up to order (u/c)^), which can be used to 
decouple the Hamiltonian constraint equation from the 
momentum constraint equations. In the ADMTT gauge 
the 3-metric is conformally flat up to order {v/c)^, at or- 
der (w/c)** deviations from conformal flatness enter. The 
extrinsic curvature up to order (v/c)^ is simply of Bowen- 
York form jlj], with correction terms of order [v/cf. 

We will use the York-Lichnerowicz conformal decom- 
position 1^ and use the PN data as the freely specifiable 
data. We numerically solve for a new conformal factor ^' 
and the usual correction to the extrinsic curvature, given 
by a vector potential . The new extrinsic curvature 
and the 3-metric multiplied by are then guaranteed to 
fulfill the constraints. The real problem in this approach 
is to find a numerical scheme which can deal with the di- 
vergences in the PN data at the center of each black hole. 
The most serious divergence occurs in the PN conformal 
factor ipp^ oi the conformally flat part of the 3-metric. 
We therefore rescale the PN data by appropriate powers 
of ippN to generate a well behaved 3-metric. If we then 
use the conformally rescaled data as the freely specifl- 
able data and make the ansatz that the new conformal 
factor 4" is the PN conformal factor tApat plus a finite cor- 
rection u, we arrive at elliptical equations which can be 
solved numerically. The splitting of the new conformal 
factor into ^' = ipp at -I- m is very similar to the puncture 
approach pq|, except that in our case the momentum 



constraint has to be solved numerically as well. 

Let us point out several issues that arise in the con- 
struction of solutions to the constraints of the full theory 
based on PN data. First of all, the accuracy of the PN ap- 
proximation increases with the separation of the binary, 
and the same is therefore true for the numerical data. 
Second, PN theory typically deals with point particles 
rather than black holes. One has to somehow introduce 
black holes into the theory, which leads to a certain arbi- 
trariness of the data near the black holes. We make the 
specific choice contained in Note that since we are 
solving elliptic equations, the data near the black holes 
affect the solution everywhere. Third, some of the PN 
expressions that we use are near zone expansions which 
are invalid far from the particles. This means we have 
data only in a limited region of space. 

Finally, the reader should be aware of the following ba- 
sic feature of the York procedure to compute initial data. 
Given valid free data, which in our case is derived from 
the PN data, the procedure projects the data onto the 
solution space of the constraints. This projection maps 
the PN data somewhere, but is the end point better than 
the starting point? We have to make sure that we do not 
loose the advantage of starting with PN data over, say, 
simply using PN orbital parameters in the conformally 
fiat data approach. After describing and resolving sev- 
eral technical issues in the construction of our data set, 
we will therefore (i) quantify the 'kick' from PN to fully 
relativistic data, and (ii) suggest a concrete method for 
improving the results of our straightforward first imple- 
mentation. 

A. Notation 

We use units where G = c = 1. Lowercase Latin in- 
dices denote the spatial components of tensors. The co- 
ordinate locations of the two particles are denoted by 
(xi,?/i,zi) and {x2,y2,Z2)- We define 

rA ■= \/{x- xaY + {y- VaY + {z- zaY- (1) 

and 

n\:= [x - XA,y - VA.z - ZA)/rA, (2) 

where the subscript A labels the particles. Furthermore 
we introduce 

ri2 ■= \J {xy - X2f + {y\ - 2/2)^ + {z\ - Z2Y (3) 

to denote the separation between the particles. All terms 
carrying a superscript TT are transverse traceless with 
respect to the flat 3-metric bij. 

II. THE PN EXPRESSIONS FOR 3-METRIC 
AND EXTRINSIC CURVATURE 

Our starting point is the expressions for the PN 3- 
metric gf^^ and the PN 3-momentum TTp^ computed in 
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the ADMTT gauge ||||. The ADMTT gauge is specified 
by demanding that the 3-metric has the form 



no longer present. If we insert Eq. into Eq. 
expand in e we obtain 



and 



4" 



and that the conjugate momentum fulfiUs 



0. 



(4) 



(5) 



We exphcitly include the formal PN expansion parame- 
ter e ~ t;/c in all PN expressions, a subscript in round 
brackets will denote the order of each term. When a PN 
term is evaluated numerically, e is set to one. 

We start with the PN expression for the 3-metric p5|| 



where the conformal factor of PN theory is given by 
^PN^l + \ (e'0(2) + £''0(4)) + 0(e')- 



(6) 



(7) 



Using the expressions for 0(2) ^-nd (/)(-4') given in p5| we 
see that the conformal factor ippN can be written in the 
simple form 



PN 



= 1 



V — 



A=l 



(8) 



where the constants Ei and E2 depend only on the 
masses toi, m2, the momenta pi, p2 and the separation 
ri2 of PN theory. They are given by 



Ea — e iriA + e 



mim2 



\2mA 2ri2 



(9) 



and can be regarded as the energy of each particle. 

Note that the PN 3-metric is singular at the location of 
each particle, since 0(2), 0(4) and ^^^(4) all go like ~ 1/rA 
as particle A is approached, and hJ^J^s^ is regular. This 
means that the strongest singularity is in V'ptv ~ 1/^1 
and that the ippN term dominates near each particle. 
Hence near each particle the 3-metric can be approxi- 
mated by 



Ea 
2rA 



1 + 77^ ^^3+0{l/ri), 



(10) 



which is just the Schwarzschild 3-metric in isotropic co- 
ordinates. For — !■ we approach the coordinate sin- 
gularity that represents the inner asymptotically flat end 
of Schwarzschild in isotropic coordinates, which is also 
called the puncture representation of Schwarzschild. This 
shows that if we write the 3-metric as in Eq. , we actu- 
ally do have a black hole centered on each particle. This 
is non-trivial since PN theory in principle only describes 
particles. 

On the other hand, if we expand the conformal factor 
in Eq. (g), the puncture singularity of Schwarzschild is 



1 + e 2*^(2) 



(4) 



32^(2) 



+e ^ij(4) 



■e%(5)+0(e^): 



which goes like 



PN 
9^J 



( const \ 



(11) 



(12) 



near each particle. One necessary condition for a black 
hole is the presence of a marginally trapped surface, and 
while the Schwarzschild metric in isotropic coordinates 
has a minimal surface at radius M/2, the term in l/r\ in 
( |l^ ) leads to a minimum in area at radius zero (ignoring 
the extrinsic curvature terms). Therefore the particle is 
not necessarily surrounded by a horizon. 

From now on we will use the 3-metric of as writ- 
ten in Eq. (||), without expanding V'pAr in in order to 
make sure that we have black holes in our data. The 
puncture coordinate singularity has replaced the point 
particle singularity. This choice is somewhat ad hoc, but 
since PN theory is not valid near the particles anyway, 
we have to make some choice, and putting in black holes 
as punctures seems natural. 



The determinant of gf^ is 



..PN 



^JPN ■ 



(13) 



since S'^-'hJ^ = 0. 

The PN expansion for the conjugate momentum is |2£ 



where 



and 



\5) --2*^(2)^3) + 2 ('^(2) '^(S)) 



^(5) = 2 m) + 2('^(2)^(3)) ■ 



(14) 



(15) 



(16) 



As in the case of the 3-metric it turns out that TTp^^ 



Eq. ( PI ) is singular, since Ti"^^^, tt^^^ and t^I^-^"^ all diverge 
at the location of each particle. But all these singularities 

H) 



in ndjy up to 0(e ) can be removed by rewriting Eq. (|l^ 

as 



'PN 



PN 



TT 



^ ^(3) + ^ 9% 



z-ij \TT 



(4) +e (0(2)^(3)) 



(17) 



which can be verified to agree with Eq. by re- 

expanding ijjpN as in Eq. (R) and keeping only terms 
up to O(e^). Hence all singularities can be absorbed by 
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4) and TT^gj are valid 

ITT 

y"(4)' y(4) 



the conformal factor, which is the basis for the puncture 
method in general |@, |l^ . 

Note that explicit expressions for 0(2) i 4'{4) ^^^d tt 
can be found in e.g. |^ or |2^. In addition Ohta et 
al also give an expression for the lapse up to 0{e^) 
and for the shift up to O(e^). The explicit expressions 
for /^^■T/^, h^'^.s, and ttH, , however, we obtained from 
Jaranowski and Schafer in a Mathematica file. 

It should also be noted that the analytic expressions 
p5[ used for the PN terms 

everywhere, while the expressions used for h?"^ 
and ^^("5) are near zone expansions 

The near zone expansion is valid only for r <C A 
■Kyjr^/{mi +TO2), where r is the distance from the par- 
ticle sources and A is the wavelength. In principle hj^^ 
should be computed from a wave equation, but in the 
near zone this equation can be simplified by replacing 
the d'Alembertian by a Laplacian. This is exactly what 
Jaranowski and Schafer |^ do to arrive at the expression 
for /i^"^ we use. In particular, the near zone expansion for 
is a spatially constant tensor field that just varies in 
time. So for the purpose of finding initial data it suffices 

to choose the initial time such that hJ^^-, vanishes. Thus 

■y(5) 

in all our numerical computations we will set hj^^^-^ = 0. 



Using the gauge condition (m we obtain 



^PN = g^^'^T^N - 0{e'). (18) 
The next task is to compute the extrinsic curvature 



(19) 



from the conjugate momentum TTpjv- With the help of 



PN 



Eqs. (|l3|) and (18), and using the expressions for tt 
in Eq. (|17|) we find that the extrinsic curvature can be 
written as 



"pat — 



+0(e«), 



3 ~y 
' ^(3) 



(20) 



such that the conformal factor ijjpM is factored out. The 
leading term in Eq. (po[) is of Bowen-York form, i.e. 



^(3) 



y- 



A=l 



Pa^a + pWa 



(21) 



outside the singularities. Moreover from Eq. ([ij) we have 



(23) 



so that K^pj^ can be considered traceless up to 0(e^ 



III. CIRCULAR ORBITS IN PN THEORY 

The PN expressions given in section |l| are valid for 
general orbits. Any particular orbit is specified by giv- 
ing the positions and momenta of the two particles. In 
this paper we want to consider quasi-circular orbits, since 
they are believed to be astrophysically most relevant. For 
a given separation r\i we therefore choose the momenta 
such that we get a circular orbit of post-2-Newtonian 
(2PN) theory. If we choose the center of mass to be at 
rest the two momenta must be opposite in sign and equal 
in magnitude. Also, for reasons of symmetry and P2 
for circular orbits must be perpendicular to the line con- 
necting the two particles. Next from the expressions for 
angular momentum and energy for circular orbits given 
by Schafer and Wex we find that the momentum 
magnitude for circular orbits is given by 



circ\2 

'pn) 



, M 



(Hi? 



+ £^(74 



' 12 



2 

(24) 



where M = mi + m2 and = mim2/M. If this for- 
mula for the momentum together with the separation is 
inserted into the expressions for 3-metric and extrinsic 
curvature in section ^, we obtain PN initial data for cir- 
cular orbits. There are, however, at least two ways how 
this can be done. One way is to always insert the mo- 
mentum (|2^ ) to the highest order known, even in terms 
which are themselves say of O(e^). One might hope to 
thereby improve the PN trajectory information in the ini- 
tial data. Another way is to consistently only keep terms 
up to a specified order, say up to O(e^). As an example 
let us look at the PN conformal factor given by Eqs. (|^) 
and (^). As one can see from Eq. (^, the momentum 
terms are already O(e^), so that if we insert Eq. ([24|), 
we generate terms of O(e^) and 0(e*), which should be 
dropped if we consistently want to keep terms only up 
to O(e^). We wiU see later that the ADM mass of the 
system is indeed sensitive to whether or not we drop such 
terms in the conformal factor. 

In order to compare with numerically computed ADM 
masses, we will also need an expression for PN total en- 
ergy of the system. For circular orbits it is given by 



Using that ^iTr^g^ = outside the singularities and the 
fact that the last two terms inside the square bracket of 
Eq. ( pO| ) are transverse (with respect to 5ij), we find 



10 T^ y" 



) = 0(6^) 



(22) 



E'l^^'^ = M- 



2^ 



A-7 

M 



9 + 20 



M 



M2 



M2 
8^ 



M 



4r 



12 



+ 0(6«). (25) 
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IV. SOLVING THE CONSTRAINTS 
A. The York Procedure 



Here Oj-' , Kl'' and ■05 are the PN expressions (^), (|2C 
and (g) with all terms of O(e^) or higher dropped. 
For (jij we choose the conformally rescaled metric 



The PN expressions for the 3-metric and the extrinsic 
curvature as given in Eqs. and (|2^) do not fulfill the 
constraint equations of general relativity. In order to find 
a 3-metric and extrinsic curvature which do fulfill the 
constraints, we now apply the York procedure to project 
the PN 3-metric and extrinsic curvature onto the solution 
manifold of general relativity. In this procedure we freely 
specify a 3-metric gij, a symmetric traceless tensor A^^ 
and a scalar K. We then solve the constraint equations 

= V^* - — ^^/ST^ 
8 12 

+ + LW^^J)(^A''i + LW''')g,kgji (26) 

8 

and 



(27) 



for and . Here V and R are the covariant derivative 
and Ricci scalar associated with the 3-metric gij , LW^^ = 
V'W^ + V^W^ - |5*^■VfeI^^ and hiM^ = V ^LW^K 
Then 



(28) 



and 



1 



i^y ^ -f LW'^) + -g'^K (29) 

3 

with g^^ being the inverse of g^ will satisfy the constraints 
of general relativity. 



- i'l^l = h + [h^m + '"^m) ' (33) 

which has the advantage of being regular near the black 
holes. We also conformally rescale the extrinsic curvature 
and pick 



- *7 

= -^3) - 



9iK,, 
(34) 



where — g^ K^^ . Finally, since we only consider terms 
up to order and because KpN = O(e^) we choose 

X = 0. (35) 

The metric is regular near the black holes. If 
denotes the distance to the singularity, we have 



and ft.^.'^^) 4" /i^Jg-, ^I/taso that 

g^j 5,j + 0{r\). 



(36) 



(37) 



This means that Christoffel symbols and Ricci scalar 
computed from the 3-metric gij go as 



and 



R^O{rA). 



(38) 



(39) 



B. Application of the York Procedure to the PN 
data 



We also have 



(40) 



The idea is to base the freely specifiable quantities gij , 
A*-' , and K on the PN 3-metric, the traceless part of the 
PN extrinsic curvature and the trace of the PN extrinsic 
curvature. The specific PN expressions we use are 



4 



f ,TT 



%-(5) 



and 



10 



~ij -.TT 



/ly(4) + (<^(2)^(3)) 



with 



1 

2~i 
1 
2^ 



9 

Pi 



mim2 



2toi 2ri2 



1712 



_P2_ 

2to2 



miTO2 
2ri2 



(30) 



(31) 



(32) 



and thus 



and 



K, ^ 0{r\) + 0{r\) 



A^' ^ 0{l/r\) + 0{l/r\). 



(41) 



(42) 



So except for A^^ and 05 all quantities are well behaved 
near the black holes. 

The remaining problem is to solve ( p6| ) and ( |2^ ) numer- 
ically. Since the PN metric is an approximate solution 
it is clear that « -05 and hence that ^ will diverge 
near the black hole, which of course is problematic when 
V^^f ^ 0{l/r^) is calculated by finite differencing in nu- 
meric computations. In order to overcome this problem 
we make the ansatz 



(43) 



6 



which in the case of the original puncture data suffices to 
regularize the constraint equations . With this ansatz 
Eq. (pq) becomes 







1 



+ ^,-7(^^.- + Ziy^^Ki^-' + Lw''')g,k9ji, (44) 



where the term 



(45) 



has been subtracted. This term vanishes analytically 
away from the punctures and it is numerically advan- 
tageous to use it to cancel the corresponding term in g"^^ . 
Using Eqs. (|6|), (^), (|3|) and (|2|) one can check that 
all terms in Eq. (H) are finite. Furthermore we split A''^ 
into the two parts 



and the York procedure would still yield a solution to the 
constraints. Each of these different starting points will in 
general yield different results for gij and Kij depending 
on Q. The solution for gij and Kij becomes independent 
of Q only if K^-' already fulfills the momentum constraint, 
which is not the case for the PN expressions. As an 
example of this freedom we expand il in e and choose 



n = l + e^Q + 0{e'^). 



(52) 



Due to the absence of O(e^) terms in we obtain the 
simple result 



0{e')fgi 



4.TT 



(53) 



and 



1:tT 



and 



(46) 



(47) 



^PN 



(l + e'*0 + O(e«))-i° 



PN 



"-PN 



0(e«)-(54) 



so that A'^ = A| 



The advantage of splitting A^^ 
in this way is that, analytically. 



We see that gfj^ and K^-' 
only in the factor 



pjY differ from gfj^ and Kp^ 







away from the punctures. Using Eq. (| 
equation (E7h simplifies to 



(48) 

the constraint 



V.Al^O. 



(49) 



Eqs. (|4J) and (^9|) now can be solved numerically for u 
and given the boundary conditions that u and 
W' for r — > oo. There are no additional boundary 
conditions at the punctures, rather we assume that there 
exists a unique solution for which u and are at 
the punctures, which has been proven to be the case for 
the simpler example considered in pM- 



Ambiguities in the application of the York 
procedure 



PN 



PN 



(55) 



This shows that an overall conformal rescaling hy U, — 
1 + e^Q can be understood as a shift (by e^Q) in the PN 
conformal factor. 

Furthermore note that any 3-metric gij and extrin- 
sic curvature ATy constructed by the method explained 
above are in general different from the PN expressions for 
3-metric and extrinsic curvature. If one assumes that the 
PN expressions are valid and thus astrophysically realis- 
tic (at least in a certain regime), one can aim to minimize 
the difference between gij and Kij and the PN expres- 
sions in this regime. We will later show that the scaling 
in Eq. ( |52| ) can be used to improve gij such that the ADM 
mass of the system after the York procedure is close to 
what is predicted by pure PN theory in the regime where 
PN theory is valid. 



Note that the York procedure explained above was ap- 
plied to the conformally rescaled quantities gij and A^^ . 
There is a priori no reason for using and A^^ . In prin- 
ciple we could have also started directly with gfj^ and 

Kpj^ or with gfj^ and Kpj^ scaled by any function fl, 
i.e. with 

gr/ = f^^^r (50) 
K'^^ = 1^-1° A7"^ (51) 



V. NUMERICS 

We now demonstrate that our method for solving the 
constraints in Eqs. (^^ and (^9|) leads to convergent nu- 
merical solutions. We use second order finite differencing 
together with a multigrid elliptic solver (BAM_Elliptic 
in Cactus [|3|). All grids have uniform resolution. The 
two black hole punctures are always staggered between 
grid points on the finest grid in the multigrid scheme. 
Since we absorb all diverging terms in the conformal fac- 
tor the solutions u and lU* of Eqs. (^ and (^) are 
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FIG. 1: Hamiltonian constraint violation for a black hole sep- 
aration of ri2 = 8M. The Hamiltonian constraint of pure 
PN data is much larger than the Hamiltonian constraint after 
solving (i.e. applying the the York procedure). We numer- 
ically solve for three different resolutions h. The inset is a 
blow up of the central region, which shows that our numeri- 
cal scheme is second order convergent as expected. 



regular everywhere, so that no black hole excision or in- 
ner boundary conditions are needed. As outer boundary 
conditions we use Robin conditions, i.e. we assume that 
u (X 1/r and oc l/r, where r is the distance to the 
center of mass. In the case of the vector potential this is 
a simplifying assumption that works reasonably well in 
practice. 

For the numerical work in this paper we consider non- 
spinning equal mass binaries with their center of mass 
at rest at the origin. The binaries are in quasi-circular 
orbits in the sense that we use Eq. ( p4| ) to set the mo- 
mentum of the two black holes before solving the con- 
straints. The two black holes are on the y-axis, such 
that their momenta point in the positive and negative x- 
directions, resulting in an angular momentum along the 
z-direction. Fig. [| shows the Hamiltonian constraint vi- 
olation of pure PN data (dashed line), i.e. before solving 
the constraints, as well as the Hamiltonian constraint 
after solving at three different resolutions h. After the 
elliptic solve the constraint equations (44) and (^) are 
satisfied to within a given tolerance of 10^"'^° in the 12- 
norm, but to study convergence we show the ADM con- 
straints computed from gij and Kij . The two black holes 
are at y = ±4. One can see that the constraint violation 
after the York procedure is much smaller than the con- 
straint violation of pure PN data. The inset in Fig. |l| is 
a blow up of the center and shows second order conver- 
gence to zero in the Hamiltonian constraint after solving. 
We also observe second order convergence to zero in the 
momentum constraint. As an example we show the y- 



FIG. 2; The momentum constraint for a separation of ri2 = 
8AI. We observe second order convergence in the resolution 
h after solving. The momentum constraint violation of pure 
PN data is larger than after solving. 




FIG. 3: The solutions of u and along the y-axis for a 
black hole separation of ri2 — 8M. For comparison we also 
show ippN, which diverges at y = ±4. 



component of the momentum constraint in Fig. ^ We 
see that pure PN data violates the constraints. In Fig. ^ 
we plot the solutions u and along the y-axis, which 
contains the black holes. As expected they are regular, 
unlike 4'pn which diverges at the black hole locations of 

As expected, after applying the York procedure g.^j and 
K^^ are different from the pure PN expressions gfj and 
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FIG. 4: Components of the 3-metric and extrinsic curvature 
for a black hole separation of ri2 = 8M. The data are shown 
before (dashed lines) and after applying the York procedure 
(solid lines). The components of the 3-metric change on the 
order of ^ 1%. 



TABLE I: Selected components of the 3-metric, extrinsic cur- 
vature and hjj^^-^ at the point x = 0, y — 12. 2M, z — for two 
black holes located on the y-axis ai y = ±5.2M. The change 
in the 3-metric induced by solving the constraints without 
first rescaling tppN has about the same magnitude as the PN 
corrections at 0(e*). The data here are computed by incon- 
sistently keeping all higher order momentum terms in ippN- 



FIG. 5: PN energy of Eq. ( p5| ) and ADM masses before 
and after solving (i.e. applying the York procedure) versus 
coordinate separation ri2 along the PN inspiral sequence. The 
data here were computed by keeping all momentum terms in 
ijjpN, without consistently dropping higher order terms. In 
this case the ADM mass of pure PN data does not agree 
well with the PN energy. The ADM mass after solving (with 
q — 0.0) increases on the order of ~ 1%, when compared 
to the ADM mass of pure PN data. Furthermore the ADM 
mass after solving increases with decreasing separation, which 
is physically not acceptable. For comparison we also show the 
ADM mass of two puncture black holes along the PN sequence 
with constant bare masses, which show a similar increase in 
ADM mass. 



PN value 
(up to Oie"^)) 


Value after 
solving {q = 0) 


relative 
difference 


g^^ = 1.21866 
K^y^ = -0.0022341 


g^^ = 1.22285 
K^y = -0.0022617 


= 0.0034 
""Ip^ = -0.012 


PN metric 
(up to 0(6^)) 


TT term in metric 
of 0(e*) 


relative size of 
0(e*) correction 


gfj^ = 1.21866 


hn(4) = 0.00443 


= 0.0036 



Fig. ^ shows a comparison of several components 
of the 3-metrics tp^^^gij and tppj^gfj^ . As one can see, 
the components of gij exhibit an increase on the order 
of ~ 1% when compared to gfj^ ■ The same conclu- 
sion is reached by looking at Tab. |, which shows the 
3-metric and extrinsic curvature before and after apply- 
ing the York procedure. Furthermore Tab. | shows that 
the increase in the 3-metric due to applying the York 
procedure has about the same order of magnitude as the 
PN corrections at O(e^). Since this happens in a region 
far enough from the particles that PN theory can ac- 
tually be trusted to give realistic values, it means that 



solving the elliptic equations introduces significant dif- 
ferences between gij and gfj^ in the outer region due 
to changes in the inner region. Before we suggest how 
this problem can be addressed, let us also consider the 
ADM mass of the system, which is a coordinate invariant 
quantity. 

We compute the ADM mass along PN inspiral se- 
quences constructed from PN circular orbits with differ- 
ent radii. Along such a sequence the bare masses mi and 
TO2 are kept constant and the momenta are computed 
from Eq. (|2j) for circular orbits. Fig. || shows the nu- 
merically computed ADM mass of pure PN initial data 
(dashed line), the ADM mass of the data obtained af- 
ter applying the York procedure (long dashed line), as 
well as the PN total energy (dotted line) of Eq. (p5|). In 
Fig. H and the following figures we plot data for ri2 be- 
tween 1 and 2QM. But note that it has to be expected 
that the PN data becomes inaccurate for small ri2, for 
example for ri2 « 4M where the black holes are close to 
the fiducial ISCO of the PN data. 

In Fig. |[ we again observe an increase of ~ 1% in the 
ADM mass after applying the York procedure. A further 
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mass/M PN conformal factor 




FIG. 6: PN energy of Eq. (ESI) and ADM masses versus coor- 
dinate separation ri2 along the PN inspiral sequence. Shown 
are the ADM masses before and after applying the York pro- 
cedure with both g = and q — 0.65. Here all data are 
computed by consistently keeping momentum terms in tppN 
only up to Newtonian order. The ADM mass of pure PN data 
now agrees better with the PN energy. The York procedure 
with q — 0.0 again increases the ADM mass on the order of 
~ 1%, when compared to the ADM mass of pure PN data. 
The ADM mass after solving with q = 0.65, however, does 
not change very much and it also closely follows the PN en- 
ergy down to ri2 ~ 6M. Furthermore until ri2 ~ 5.6M it is 
physically reasonable since it decreases with decreasing sepa- 
ration. For comparison we also show the ADM mass curve of 
rescaled PN data (with q = 0.65). These data, however, have 
no direct physical significance. 



problem is that none of the numerically determined ADM 
masses in Fig. |^ agrees very well with the PN energy (25 ). 
This problem stems from the fact that the PN initial data 
in Fig. 1^ have been obtained by inserting the momentum 
( p^ ) as it is into the expressions for 3-metric and extrinsic 
curvature of Sec. |l| without consistently dropping terms 
of O(e^) or higher. Since all PN corrections to the mo- 
mentum are positive, the main effect of this inconsistency 
is to increase ipPN given by Eqs. (||) and (||). The result is 
that the numerically computed ADM masses before and 
after applying the York procedure show physically unac- 
ceptable behavior: (i) the ADM mass of pure PN data 
approaches the PN energy (j2^) only very slowly at large 
separations, and (ii) the ADM mass of the data after ap- 
plying the York procedure monotonically increases with 
decreasing separation. This is physically not reasonable 
because the system is supposed to loose energy due to 
the emission of gravitational radiation. For reference the 
ADM mass (dot dashed line) for a sequence of two black 
hole punctures with constant bare masses and with the 
same PN momentum (|2^) is also shown in Fig. ||. Along 



FIG. 7: The conformal factors TppN and i/ipjv ~ ipPN + e^Q, 
before and after rescaling with q — 0.65 for ri2 = 8M. The 
difference between tpPN and tppN is small. 



this sequence the ADM mass of the punctures also un- 
physically rises with decreasing separation, which is not 
surprising since the assumption of constant bare masses 
for punctures ignores the growing contribution of u to 
the conformal factor with decreasing separation of the 
punctures. In all cases studied by us the solution u of 
Eq. ( p[ ) is indeed positive, which translates directly into 
an increase in the mass. 

Of course, the question is how we can improve our data 
so that its behavior is physically more realistic. One can 
argue that part of the additional energy is tied to an 
increased local mass of the individual black holes. In 
fact, for constant bare masses there is a strong growth 
in the apparent horizon masses. A standard approach is 
therefore to rescale the bare masses to keep the apparent 
horizon mass fixed and to define a binding energy by 
subtracting the apparent horizon masses from the total 
mass, e.g. We plan to compute apparent horizon 

masses for a future publication. However, in general it 
is not possible to unambiguously define a local mass for 
general relativistic data, and the accuracy and validity 
of the estimate for the binding energy therefore depends 
on, for example, how close the black holes are. 

As an alternative we have experimented here with a 
mass correction that is tied to properties of the PN ap- 
proximation. As a first step let us keep momentum terms 
of Eq. ( |24| ) in the PN conformal factor ipPN (see Eqs. (^) 
and (^)) only up to the appropriate order and to consis- 
tently drop all terms of O(e^) and higher. This amounts 
to just using the first Newtonian term of the momen- 
tum ( p^ ) in t/jpN- The results are shown in Fig. ^. The 
ADM mass of pure PN data (dashed line) now much bet- 
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ter approaches the PN energy for large separations. Yet, 
the ADM mass after simply applying the York proce- 
dure (long dashed line) still shows an increase of order 
~ 1% when compared to pure PN data. If we want more 
physical mass curves we have to prevent this increase 
by preventing the increase in the conformal factor. We 
will take advantag e of th e freedom in the York procedure 
mentioned in Sec. IV C and use the conformal rescaling 
of E g. (|5^ ) before applying the York procedure. From 
Eq. ( pq ) we see that then the overall conformal factor 
becomes 



mass/M 



v]/ = -0 



PN 



PN 



(56) 



Hence, if we choose an appropriate Q, we have a chance 
of compensating u such that ^ « '^Ppn at least in the 
region far from the black holes where PN theory is valid. 

Now, in the limit of ri2 oo the pure PN data we 
use as a starting point represent two Schwarzschild black 
holes at rest (in isotropic coordinates). Thus u is zero 
for infinite separation and we therefore expect that u 
goes like u oc l/ri2 for large ri2. On the other hand 
we also have u cx l/r so that we expect that u is well 
approximated by 



N 
ri2r 



(57) 



for large r, where N is some numerical constant. Since 
we want Q + m ~ 0, we choose 



2ri2 V2ri 



1 



1 
2^ 



(58) 



Here g is a free parameter, which has to be chosen such 
that Q + u ~ for large separations. 

It turns out that for q = 0.65 we get physically more 
reasonable mass curves in the regime where PN theory is 
expected to be valid. The solid line in Fig. ^ shows the 
ADM mass obtained for different separations if we apply 
the following extended York procedure: (i) start with the 
pure PN initial data, (ii) rescale i/'PAf using Eqs. ( ^5|) and 
( p8| ) with q — 0.65 and (iii) apply the standard York pro- 
cedure to the rescaled quantities. As we can see the ADM 
mass (solid line) closely follows the PN energy (dotted 
line) in the region where we expect PN theory to be valid. 
Furthermore for separations greater than ri2 ~ 5.6M 
the ADM mass decreases with decreasing separation as 
it should. For smaller separations the ADM mass again 
increases. In the literature this minimum has often been 
interpreted as the location of the innermost stable circu- 
lar orbit (ISCO). Note, however, that the PN expressions 
which we used up to O(e^) are probably close to break- 
ing down around Mjrxi = 1/5.6 ~ 0.2, so that the ISCO 
location may not be very accurate. Also the location of 
the minimum can be shifted if we use higher order terms 
in the rescaling of ippN, i-e. if we use 
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FIG. 8: PN energy of Eq. (^, ADM mass of pure PN 
data and ADM mass after solving (with q = 0.65) versus the 
PN angular velocity (60|). The PN energy has a minimum 
near Mlopn ~ 0.1, which is often interpreted as the ISCO. 
We see that the ADM mass after solving (with q — 0.65) 
closely follows the PN energy until Mujpm ~ 0.05. Then near 
MujpM « 0.06 it has a minimum which could be regarded 
as the ISCO. One has to keep in mind, however, that the 
ambiguities in the York procedure in principle allow us to 
shift the location of this minimum. 



The extra Q' term will have no influence in the limit of 
large distances, but it will influence the mass curves at 
small separation and thus we can move the minimum. 
Again one could introduce a one-parameter family of Q' 
terms and flt the parameter such that the ADM mass 
curve has the minimum at the same place where the 
PN energy (|2^) has a minimum. We decided not to do 
this since the PN energy itself may not be very reliable 
near its minimum. For comparison. Fig. ^ also shows 
the ADM mass curve (dot dashed line) for the PN data 
rescaled by Q with q = 0.65, but without applying the 
York procedure. This curve has no direct physical mean- 
ing, but we can see that it can be obtained from the 
curve for pure PN data (dashed line) by a downwards 
shift. Fig. ^ shows the PN conformal factor before and 
after rescaling with q — 0.65. We see that the change in 
'0pAr is rather small. 

All the masses so far are plotted versus the coordinate 
separation ri2. Fig. || shows the PN energy (dotted line), 
the ADM mass of pure PN data (dashed line) and the 
ADM mass of data obtained after rescaling with q = 0.65 
and applying the York procedure (solid line), versus the 
PN angular velocity ujpNi computed for circular orbits 
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TABLE II: Selected components of the 3-metric, extrinsic cur- 
vature and ^^"(4) at the point x = 0, y = V2.2M, z — Q for two 
black holes located on the y-axis at ?/ = ±5.2M. The change 
in the 3-metric induced by solving the constraints after first 
rescaling -i/ipjv (with q — 0.65) is much smaller than the the 
PN corrections at O(e^). The change in the extrinsic curva- 
ture due to solving, however, does not depend much on q and 
is about the same whether or not we use the rescaling with 
q = 0.65. Here we have included only Newtonian momentum 
terms in ■i/'pjv, in order to have a consistent expansion in e. 



PN value 
(up to 0{e^)) 


Value after 
solving (g — 0.65) 


relative 
difference 


g^^ = 1.21738 
K^y'^ = -0.0022353 


g^^ = 1.21783 
K^y = -0.0022673 


= 0.00037 
-^^mr^ = -0.014 


PN metric 
(up to 0{e^)) 


TT term in metric 
of O(e^) 


relative size of 
0(6*) correction 


g^^ = 1.21738 


/i^4) = 0.00443 


= 0.00364 



from 



64(ri2/M)3 
(1 + 2ri2/M)6 ^ 

5 /i fi^ 



M 



+ ^ - ■ (60) 



M 

M' 
ri2 



Note that upN in Eq. (60) is written such that wp^r is 
exact up to all PN orders in the limit of /i/Af — + 0. For 
^/M > Eq. ( |60| ) is accurate up to 2PN order. It should 
be kept in mind, however, that lopn probably is not ex- 
actly equal to the true angular velocity after applying 
the York procedure. Yet our numerical approach does 
not immediately yield an angular velocity which could 
be used in place oftupN- 

From Fig. || we see that the approximate ISCO of PN 
theory computed from the 2PN energy is near MuopM — 
0.1, while the ISCO minimum of our data (after apply- 
ing the extended York procedure with q — 0.65) is near 
MujpM = 0.06, which is very close to the ISCO of test 
particles in Schwarzschild. Also note that the ADM mass 
of pure PN data (dashed line) does not have a minimum 
at all. 

In Tab. |l| we compare some components of the 3- 
metric and extrinsic curvature of pure PN data with the 
corresponding quantities obtained after rescaling with 
q — 0.65 and applying the York procedure. The change 
in the 3-metric induced by solving the constraints after 
first correcting -0pjv (with q = 0.65) now is much smaller 
than the the PN corrections at 0{e'^). The change in the 
extrinsic curvature due to solving, however, is nearly the 
same whether or not we use the rescaling with q = 0.65. 

The question arises if the solutions gij and K^^ with 
q — 0.65 are astrophysically more realistic then the pure 
PN solutions gfj^ and Kpj^. We argue that this is indeed 

the case since gij and K^^ with q — 0.65 are close to g[j^ 



and Kpj^ in the far region where PN is accurate, 



but 



in addition do fulfill the constraint equations of general 
relativity. Furthermore the ADM mass curve for gij and 
K^^ with q = 0.65 is closer to the PN energy ( p5| ) than 
the ADM mass curve of the pure PN solutions gfj^ and 



VI. DISCUSSION 

For the first time, we have derived fully relativistic 
black hole initial data for numerical relativity, starting 
from 2PN expressions of the 3-metric and extrinsic cur- 
vature in the ADMTT gauge. We have used the York 
procedure, and any procedure for projecting the PN data 
onto the solution manifold of general relativity will intro- 
duce changes to the PN data. The larger the violation of 
the constraints by the PN data, the larger the change in 
the solution process will be. In principle one may loose 
the PN characteristics that distinguished the PN data 
from other approaches in the first place. 

As we have seen in Sec. the size of these changes de- 
pends on how exactly we employ the York procedure for 
the projection. We find that the extended York proce- 
dure (with q — 0.65) yields acceptably small changes, so 
that if the PN data we started with are astrophysically 
realistic, the data after solving the constraints should 
still be astrophysically relevant. In particular, our new 
PN initial data have the nice property that the 3-metric 
and extrinsic curvature approach the corresponding 2PN 
expressions in the region where PN theory is valid, pro- 
viding a natural link to the early inspiral phase of the 
binary system. Furthermore, our approach leads to an 
easy numerical implementation with a generalized punc- 
ture method. 

We consider this work as a first step towards the con- 
struction of astrophysical initial data based on the PN 
approximation. Although we are able to remove some of 
the inherent ambiguity of the method, several directions 
should be explored. Since the PN formalism is unable to 
unambiguously provide the full information in the black 
hole region, one should examine different ways to intro- 
duce black holes. Furthermore it would seem natural 
to follow the conformal thin sandwich approach in or- 
der to obtain data that corresponds more closely to a 
quasi-equilibrium configuration, although in principle we 
rather want data for the appropriate PN inspiral rate 
than for exactly circular orbits. Note that after the solu- 
tion process it is not known how well the orbital param- 
eters correspond to quasi-circular orbits. One could use, 
for example, the effective potential method ||l5| with the 
new PN based data to determine quasi-circular orbits of 
the two black holes. 

Another direction of research is to improve the PN 
input to our method. Even though we can solve the con- 
straints for rather small separations of the black holes, 
we cannot trust the numerical data for arbitrarily small 
coordinate separation, because this is where the PN data 
we start with is probably unreliable. We have started 
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with a traditional PN approach ||25|], but there has been 
significant progress in extending the vahdity of the PN 
approximation to smaller separations through resumma- 
tion techniques ^ |36|. It is an important issue to 
study how large an intermediate binary black hole regime 
might be, where the PN approximation has broken down 
but the separation is still significantly larger than the 
separation for an approximate ISCO ||37| . 

In addition, we want to work with higher order PN 
approximations. The explicit regularization for 3PN of 
p6[ could be used as a starting point. However, our 
procedure may have to be modified because of changes 
in the conformal factor ipPN- Finally, Jaranowski and 
Schafer |^ have recently provided us with an expression 
which includes spin terms at order {v/c)'^ in the PN ex- 
trinsic curvature. In future work we intend to use these 
terms to add spin to the black holes. 

Recall that we have concentrated on the near zone. 
We plan to replace the near zone expansion of ^^(4) 
with a globally valid expression. This could be achieved 
by solving the wave equation determining hjj^^^ (see 
e.g. p9| ) numerically, without any near zone approxima- 
tions, which would be natural in a method that resorts to 
numerics anyway. If the PN inspiral trajectory is used in 
this calculation, the initial slice of our spacetime will al- 
ready contain realistic gravitational waves, with the cor- 
rect PN phasing. When this spacetime is then evolved 
numerically we might eventually be able to compute nu- 
merical wave forms which continuously match PN wave 
forms. 

This brings us to the final goal of our initial data con- 
struction, namely to use it as the starting point for nu- 
merical evolutions. As we pointed out in the introduc- 



tion, there are now numerical evolution methods with 
which we can begin to explore the physical content of any 
initial data set by evolution and by extraction of physical 
quantities such as detailed wave forms or total radiated 
energies 0, As mentioned in |^, the Lazarus ap- 
proach provides an effective method for cross-checking 
the validity of the results by choosing different transition 
times along the binary orbit in the region where a far 
limit approximation (such as the PN method) and full 
numerical relativity overlap. Only by extending the abil- 
ity of full numerical codes to accurately compute several 
orbits, will we be able to arrive at a definitive conclusion 
about the merit of different initial data sets. 
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